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An  array  of  sensors  is  receiving  radiation  from  a source  of 
interest.  The  source  and  the  array  are  in  a closed  three  dimensional 
waveguide.  The  maximum-liklihood  estimators  of  the  coordinates  of 
the  source  are  analysed  under  the  assumptions  that  the  noise  field 
is  Gaussian,  and  the  source  lies  on  a finite  linear  manifold  spanned 
by  a set  of  orthogonal  functions.  These  functions  are  the  eigen- 
functions of  the  wave  operator  given  the  boundary  conditions  defining 
the  waveguide,  and  are  called  normal  modes.  The  Cramer-Rao  lower 
bound  is  of  the  order  of  the  number  of  modes  which  define  the  source 
excitation  function. 
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Maximum  Likelihood  Estimation  of  the  Position  of  a 
Radiating  Source  in  a Closed  Waveguide 

Melvin  J.  Hinich 

Introduction 

Arrays  of  Interconnected  receivers  are  used  to  improve  the  signal- 
to-nolse  ratio  of  coherent  radiation  from  a source  located  in  the 
same  waveguide  as  the  array.  A waveguide  is  a physical  science  term 
for  the  classical  wave  operator  and  its  associated  boundary  conditions. 
In  most  examples,  the  phases  of  the  received  signals  are  adjusted  by 
digital  or  analog  methods  in  order  to  concentrate  the  received  energy 
in  a narrow  beam  with  respect  to  the  coordinate  system  of  the  array. 
This  procedure  (called  beamforming)  Is  used  to  estimate  the  direction 
of  a point  source  which  Is  radiating  coherent  energy  at  a single 
frequency  or  In  a band  of  frequencies.  Beamforming  Is  widely  used 
for  processing  arrays  In  radio  astronomy,  underwater  acoustics,  phased 
array  radars,  seismology,  and  atmospheric  physics.  The  statistical 
properties  of  beamforming  has  been  studied  by  Levin  [11],  Clay,  Hinich, 
and  Shaman  [6],  and  Green,  Kelley,  and  Levin  [7]. 

Beamforming  is  an  essential  element  of  the  maximum  likelihood 
estimator  of  source  direction  If  the  received  radiation  is  a plane 
or  cylindrical  wave.  Beamforming  will  result  In  a biased  or  very 
imprecise  estimate  of  direction  If  the  received  wave  fronts  are 
significantly  effected  by  reflections  from  the  boundaries  of  the 
waveguide,  or  are  convoluted  by  refraction  and  dispersion  In  the 
medium  of  propagation  (see  Clay  [4]).  Such  Is  the  case  when  the 
problem  Is  estimating  the  depth  of  a source  In  an  Infinite 


1 


2 


stratified  horizontal  waveguide.  Clay  (5]  rcvic’ws  signal  processing 
theory  for  such  a waveguide  using  both  horizontal  and  vertical  arrays, 
and  relates  the  statistical  models  to  physical  oceanography.  The 
maximum  likelihood  estimator  of  source  direction  using  a horizontal 
array  is  given  by  Capon  et.  al  [3],  and  Hlnlch  and  Shaman  [8].  The 
maximum  likelihood  estimator  of  source  depth  using  a vertical  array 
is  given  by  Hlnich  [9],  [10]. 

This  paper  considers  the  problem  of  estimating  the  coordinates  of 
a source  in  a closed  three  dimensional  waveguide  when  the  source 
excites  a finite  number  of  eigenfunctions  of  the  wave  operator.  The 
closure  of  the  environment  results  in  a fundamental  limitation  of 
the  precision  the  estimator  Irrespective  of  the  number  of  receivers 
used  in  the  array.  This  limitation  is  due  to  the  limited  dimensionality 
of  the  model  which  is  Induced  by  the  physics  of  the  system. 

1.  Normal  Mode  Solution  of  the  Wave  Equation 

The  steady-state  solution  of  wave  equation  for  an  infinite  hori- 
zontal waveguide  is  a special  case  of  the  general  solution  of  the 
inhomogeneous  Sturm-Liouville  partial  differential  operator  on  the 
space  12(5)  of  square  Integrable  functions  on  xcS,  a smoothly  closed 
and  bounded  subset  of  Euclidean  space  (Chapter  5,  Vladimirov  [14]). 

For  expository  purposes  we  first  limit  attention  to  one  space  dimension, 
and  set  S - [0,1].  For  x in  the  unit  Interval  [0,1],  let  f(x)  denote 
the  twice  continuously  differentiable  solution  to  the  inhomogeneous 


Sturm-Liouville  problem  Lf  ••  s(x-x^),  where  the  operator  is  defined 
by 


Lf 


3 3 , 

■ 3x  P 3x  ^ 


+ qf. 


(1) 


s(x-x^)  is  the  intensity  of  a source  at  position  p is  a given 

continuously  differentiable  positive  function  of  xcS,  q is  non- 
negative and  continuous  in  [0,1],  and  f satisfies  the  following 
boundary  conditions  given  ^ 0,  i “ 1,...,4: 


c^f(O)  - f(0)  - 0 

C3f(l)  ‘'a  “ °- 


(2) 


For  steady-state  processes  f(x)  is  the  solution  to  a general  wave 

equation  (or  a diffusion  equation)  with  the  given  boundary  conditions 

with  normalized  space  units.  The  time  component  in  the  solution  is 

separable  from  the  space  component . - 

The  Sturm-Liouville  operator  L is  self-adjoint,  its  eigenvalues 

i >2  i discrete  and  non-negative,  and  its  eigenfunctions 

^1*  ^2*  ’ ' ' complete  in  As  an  example,  consider  the 

simplest  model  of  the  ocean  as  an  acoustic  waveguide,  the  homogeneous 

compressible  fluid  waveguide  with  a rigid  bottom  and  a free  surface. 

In  the  absence  of  gravity  effects,  the  eigenfunctions  for  this 

waveguide  are  f > /2  cos  y x,  and  the  eigenvalues  are  y 
mm  m 

■ (m-f  l/2)x,  m*  1,  2,...,  where  the  units  have  been  normalized  so 
that  the  depth  of  the  guide  is  one  unit. 
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In  the  general  model,  assume  that  ^ 0 (L  is  non-singular). 


Then 


f(x) 


E 

m“l 


(s,  4>_)y_  4'_(x) 

m m m 


(3) 


where 


(s,  (J)  ) = s(x)4>  (x)dx 

m ' u m 


is  the  inner  product  between  s and  4>  . Moreover  the  eigenfunctions 

m 

are  orthonormal,  i.e.  (41  , 4>  ) ■=  1 and  (4>  » <(’  “ 0 for  m / m' . 

mm  mm 

When  L is  the  wave  operator,  the  4>  are  called  the  normal  modes. 

m 

Suppose  that  the  source  excites  only  the  first  M eigenfunctions, 

I.e.  the  source  function  lies  on  the  linear  manifold  spanned  by 

, . . .,  ^ . The  source  function  then  can  be  written 
1.  n 

M 

s(x-x  ) - } e (x  )4>  (x),  (4) 

o £ I m o D 

m**l 

where  0, , . . .,0„  are  unknown  weights  which  are  twice  continuously 
1 n 

differentiable  functions  of  the  source  position  x . For  any  scl_(S), 
the  function  can  be  approximated  to  any  preassigned  tolerance  by 
choosing  M sufficiently  large.  As  long  as  M is  finite,  it  is  clear  from 
(4)  that  the  source's  energy  does  not  emanate  from  a single  point,  i.e. 
in  physics  terminology  the  source  is  not  a point  source.  By  the  ortho- 
gonality of  the  eigenfunctions,  it  follows  from  (3)  and  (4)  that 

M 

f (x|x  ) - V 

O / * n O IB  ID 
m-1 

Thus  the  received  signal  also  lies  on  the  linear  manifold. 

' 1 n 
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Expressions  (4)  and  (5)  holds  with  ^cS  for  ScR  > and  a source 

which  excites  only  the  first  M normal  modes.  The  eigenfunctions 

satisfy  the  boundary  conditions 
m 

a(x)f(x)  + 6(x)  = 0,  (6) 

^ ^ dn 

where  jn  is  the  external  normal  vector  to  the  (piecewise  smooth)  boundary 
of  S,  and  a,  E ^ 0 are  continuous  on  the  boundary  with  a + 6 > 0. 
Moreover  the  operator  is 

L - - dlv(pV)  + q,  (7) 

where  V is  the  gradient  operator  and  dlv  is  the  divergence.  Estimating 

the  source  position  x Involves  a linear  statistical  model  whose 

*^o 

structure  is  determined  by  the  physics  of  the  waveguide.  The  eigen- 
functions of  the  operator  L constitute  a "natural"  basis  for  the  linear 
design. 

2.  A Statistical  Model 

Assume  that  the  source  and  the  receiving  array  is  in  a closed  three 
dimensional  waveguide  as  modeled  by  expressions  (6)  and  (7).  The 
finiteness  of  the  waveguide  shapes  the  structure  of  the  noise  field 
If  the  source  excites  only  the  first  M normal  modes. 

Let  denote  the  positions  of  the  receivers  In  the  array, 

and  as  before  let  denote  the  position  of  the  source  of  Interest. 

The  coordinates  (*^2*  *o2’  *o3^  ^o  unknown  parameters  which 

are  to  be  estimated  from  a sample  of  the  signals  received  by  the  array. 

Suppose  that  the  signal  has  phase  coherent  energy  In  a narrow  band 
about  frequency  f^.  Let  y(^^)  denote  the  output  from  the  1th  receiver 
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in  the  array  which  has  been  filtered  in  a narrow  band  about  For  a 

discussion  of  filtering  In  mathematical  statistical  terminology,  see 
Anderson  [1]  and  Brill Inger  [2].  If  the  signal  has  energy  in  many 
frequency  bands,  the  results  in  this  paper  apply  to  each  narrowband 
component  of  the  received  signal. 

By  the  linearity  of  the  wave  operator. 


where  e(x. ) is  the  filtered  signal  received  at  x.  from  the  ensemble  of 
noise  sources  in  S.  Assuming  that  the  phases  of  the  noise  field  in 
the  narrow  band  around  m are  Incoherent,  the  signal-to-nolse  ratio 
is  increased  by  filtering  the  received  signal  in  the  f^  band. 

Assume  that  each  noise  source  excites  at  most  N modes,  where 
N > > M.  It  then  follows  from  (4)  and  (5)  that  lies  on  the 


linear  manifold  spanned  by 


k Let  c replace  0 in 
m m 


the  linear  combination  relating  c(x)  to  the  ^ (x) . 

OB 

In  order  to  simplify  the  stochastic  structure  of  the  noise  field 
c(x),  assume  that 


N 


(8) 


where  c, , . . . , are  uncorrelated  Gaussian  random  variables  with 
1 N 


zero  mean  and  variances  Ec 


of  the  noise  observed  at  n locations  x, , 
matrix 


Thus  the  variance-covariance  matrix 
. . , ^ Is  the  n X n 


-1  -1  T 


(9) 
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where  ♦ is  the  n x N matrix  whose  (1,  mth)  element  is  4i  (x^),  F,, 

/V,  m '*'1 

is  the  N X N diagonal  matrix  with  elements  y , and  is  the  N x N 

m N 

2 

diagonal  matrix  with  elements  m (m  = 1,  . . . , N) . 

03 

The  matrix  I is  positive  seml-def Inlte  with  rank  at  most 
^x  ^ 

N if  n ^ N.  As  can  be  seen  from  expression  (8),  there  are  only  N 
independent  realizations  from  the  stochastic  noise  process  during  t!ie 
sampling  period  regardless  of  the  number  of  sensors  in  the  array. 
Consequently  there  is  no  reason  to  have  n greater  than  N.  In  practice, 
the  number  of  modes  needed  to  provide  a good  approximation  to  the 
noise  field  must  be  experimentally  determined,  as  is  the  case  for  the 
source  radiation. 

Setting  the  number  of  sensors  equal  to  N,  suppose  that  the 
sensor  locations  are  chosen  so  that  is  non-singular.  Thus, 

tT^  = ra  (lo) 

The  Fisher  information  is  the  positive  definite  3x3  matrix 


1 (x  ) - (Vf  ^(Vf  ) 
^X  /wx  -*'X 


where  Bf  is  the  N x 3 matrix  whose  1th  row  In  the  gradient 

(I f I*  ) * I f I X ) , I f (x . 1 x ) ) at  X . Since 

3x  - A.l'^o  3x  - ' 3x  - «-l 

ol  o2  o3 

M 

f(x|x  ) “ y e (x  )y  (x),  then  from  (10) and  (11)  we  have  for  all 
/ . m o m m 

N > M, 


I (x  ) - I^8(x  )]V;^lve(x^)] 


2 2 

where  (2^  is  the  diagonal  matrix  with  elements  oi^,  . . • i 

V0(x  ) is  the  M X 3 matrix  whose  m,jth  element  is  0 (x  ). 

/V  .vO  , aX  m 

oj 
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It  lb  clear  from  (12)  that  the  elements  of  I are  Independent 

/V  X • 

of  the  array  size  as  long  as  n ^ M.  This  means  that  only  M sensors 

are  needed  for  the  maximum  likelihood  estimation  of  x , provided  that 
<vo 

the  X. 's  are  chosen  to  make  the  rows  of  ♦ linearly  Independent  vectors. 

y**  i ^ 

This  result  still  holds  when  the  variance-covariance  matrix  of  the 
is  not  diagonal,  i.e.  the  are  correlated . In  this  case  the 
matrix  in  (12)  is  the  M x M variance-covariance  of  (e^^ ' 

For  each  j “ 1,  2,  3 and  m = 1,  . . . , M,  let  denote  the 

3 T 

maximum  of  It 0 (x  ) | for  x eS.  Given  a vector  a = (a,,  a_ , a.) 

' 3x  , m -^.o  ' ~o  1 2 J 

oj 

of  unit  length,  it  follows  from  (11)  by  the  Schwarz  Inequality  that 


M 

T - ? 2 2 2 

a I (x  )a  < I 0)  (B  , + B , + B _) 
/V  X /V  — , m ml  m2  m3 

m«l 


-1  2 

< 3A  B M, 


(13) 


2 2 

where  B « max(B  and  A = mln{u)  }.  Thus  the  variance  of  an  unbiased 
m.j  m — 

estimator  of  x , Is  bounded  below  by  A(3B  M)  . For  a non-diagonal 
oj 

A Is  replaced  by  the  minimum  eigenvalue  of 


3.  Maximum  Likelihood  Estimation  of 

Assume  that  the  M sensor  array  is  constructed  so  that  the  M x M 
matrix  ♦ • non-singular.  Then  the  model  can  be  written 

y “ ♦ r ^(0(x  ) + e),  (14) 

where  ^ - (y(x^) 

T 

and  e ■ (e^j^) e(Xj^))  . Since  the  e^^)  are  normally  distributed 

with  variance-covariance  matrix then  the  first  order  conditions  for 
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the  maximum  likelihood  estimator  of  x are 


(r*  - 0(x  ))'^  n"^[ve(x  )]  = o, 

^ ^ o M ^o 


where  V0(x  ) is  the  M x 3 matrix  (r 0 (x  )).  The  maximum  likelihood 

ArO  OX  , m 

A Oj 

estimator  x is  ona  of  the  solutions  to  equation  (15)  (see  Chapter  5, 

•~o 

Rao  [12]). 

It  has  been  assumed  that  0(x  ) is  twice  continuously  differentiable. 

<*»o 

g3 

In  addition,  assume  that  I- r r 0(x  )|  has  a finite  upper  bound  for 

oi  oj  ok 

all  x^eS  and  1,  j,  k = 1,  2,  or  3.  Given  the  normality  of  the  errors 

k ^ 

in  the  model,  it  follows  that  as  M -+■«>,  M (x  - x ) is  asymptotically 
normally  distributed  N(0,  I ^(x  )). 

^X 

A . Nonstationar ity  of  the  Waveguide 

It  is  often  the  case  that  the  medium  or  the  boundaries  of  the  wave- 
guide are  slightly  nonstationary  during  the  period  when  the  received 
signals  are  being  filtered.  In  order  to  model  the  effect  of  this  non- 
stationar ity  , assume  that  the  filtered  signal  received  at  x from  a 

source  at  x is  given  by  the  expression 

M 

f(x.|x  )'-  y 6„(x„)It_^'I'_(x.)  +u  (x  )],  (16) 

A/  1 W O / j II|A>0  IS  IDa^I  D 1 

2 2 

where  Eu(x. ) - 0,  Eu  (x, ) « o , and  u(x. ) Is  Independent  of  e(x  ) 

ni  Ao  X ID  A*  X 

for  each  i *•  1,  . . .,  M.  In  other  words,  expression  (5)  is  modified 
by  adding  a stochastic  perturbation  term  u(^^)  to  (see  Chapter  6, 

Tolstoy  and  Clay  [13]).  The  results  given  in  the  previous  sections  still 
hold  as  long  as  the  0 (x  ) have  the  previously  assumed  smoothness, 

ID 

but  is  replaced  by_^flj^  + where  1^  is  the  variance-covariance 
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matrix  of  the  perturbation  terms.  Assume  that  1 Is  non-singular  and 

u 

is  Independent  of  the  array  geometry.  The  perturbation  effect  raises 

the  asymptotic  variances  of  the  maximum  likelihood  estimators  of  the 

X since 
oj 


T 

a 


[V0(x  )] 


[76(x  )]a  < a”^[Ve(x  )]'^(0 

O ^ ^ /«v  O >v 


M 


+ 1 )"^[76(x  )]a 

/V-  U O 


(17) 


for  a 0 and  ^ 0.  This  result  Is  hardly  surprising.  It  Is  also 

2 

clear  that  if  is  very  large  for  large  m,  the  precision  of  the 
maximum  likelihood  estimator  is  limited  as  if  M was  bounded. 
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